Zauważmy, że zastosowane metody Eulera sprowadza się do wykonania sekwencji izogeometrycznych L2 projekcji, podobnie jak w przypadku obliczania projekcji bitmapy. Tym razem jednak naszą "bitmapą" jest stan ukłądu w poprzedniej chwili, plus zmiany wywołane przez "fizykę" modelowanego zjawiska w czasie kroku czasowego, plus zmiany wywołane siłą działającą na nasz system w czasie kroku czasowego.
Sposób traktowania równania niestacjonarnego jako sekwencje izogeometrycznych projekcji został opisany w pracy "IGA-ADS: Isogeometric analysis FEM using ADS solver" [1].
Żeby opracować solwer wykorzystujący metodę Eulera w izogeometrycznej metodzie elementów skończonych, musimy przetransformować sformułowanie silne w sformułowanie słabe.
Przemnażamy więc nasze równanie przez funkcje testujące
\( (u_{t+1},w) = (u_t + dt * \mathcal{L}(u_t)+ dt* f_t,w) \)
Do aproksymacji stanu naszego systemu w danej chwili czasowej użyjemy kombinacji liniowej funkcji B-spline. W tym celu wybieramy bazę dwuwymiarowych funkcji B-spline, określając wektory węzłów w kierunku osi układu współrzędnych. Dla ustalenia uwagi możemy wybrać dwuwymiarową bazę funkcji B-spline drugiego stopnia
\( \{B^x_{i,2}(x)B^y_{j,2}(y)\}_{i=1,...,N_x;j=1,...,N_y} \)
Będą one stosowane do aproksymacja symulowanego pola skalarnego w aktualnej chwili czasowej
\( u(x,y;t+1)\approx \sum_{i=1,...,N_x;j=1,...,N_y} u^{t+1}_{i,j } B^x_i(x)B^y_j(y) \)
Podobnie do testowania użyjemy funkcji bazowych B-spline:
\( \{ B^x_k(x)B^y_l(y) \}_{k=1,...,N_x;l=1,...,N_y } \)
Nasze równanie wygląda teraz identycznie jak problem dwuwymiarowek izgeometrycznej L2 projekcji
\( \sum_{i=1,...,N_x;j=1,...,N_y} u^{t+1}_{i,j } (B^x_i(x)B^y_j(y),B^x_k(x)B^y_l(y))=RHS \quad \forall k=1,...,N_x; l=1,...,N_y \)
Nasza prawa strona nie jest jednak bitmapą, ale projekcją sumy trzech elementów:

1. Stanu naszego układu w poprzedniej chwili czasowej \( (u_t,w)= \sum_{i=1,...,N_x;j=1,...,N_y} u^t_{i,j } (B^x_i(x)B^y_j(y),B^x_k(x)B^y_l(y)) \) (również przemnożonego przez funkcje testującą i odcałkowanego po obszarze). Oczywiście do przedstawiania stanu naszego układu w poprzednim kroku czasowym również wykorzystujemy kombinacje liniową funkcji bazowych B-spline \( u(x,y;t)=u_t=\sum_{i=1,...,N_x;j=1,...,N_y} u^t_{i,j } B^x_i(x)B^y_j(y) \)
2. Zmian wywołanych przez "fizykę" symulowanego zjawiska, w czasie kroku czasowego. Zmiany te policzone są poprzez zastosowanie operatora różniczkowego reprezentującego modelowane zjawisko do stanu układu w chwili poprzedniej. Stan układu oczywiście reprezentowany jest przez kombinacje liniową funkcji B-spline. Tak więc operator różniczkowy stosowany jest do funkcji B-spline \( (dt * \mathcal{L}(u_t),w)=\sum_{i=1,...,N_x;j=1,...,N_y} dt*u^{t}_{i,j } ({\cal L} (B^x_i(x)B^y_j(y)),B^x_k(x)B^y_l(y)) \). Na przykład dla problemu transportu ciepła mamy \( \sum_{i=1,...,N_x;j=1,...,N_y} dt*u^{t}_{i,j } (\frac{\partial^2 (B^x_i(x)B^y_j(y))}{\partial x^2}+\frac{\partial^2 (B^x_i(x)B^y_j(y))}{\partial y^2},B^x_k(x)B^y_l(y)) \). Dzięki sformułowaniu słabemu możemy odcałkować przez części \( \sum_{i=1,...,N_x;j=1,...,N_y} dt*u^{t}_{i,j } \left( (\frac{\partial (B^x_i(x)B^y_j(y))}{\partial x}B^y_j(y),\frac{\partial (B^x_k(x)B^y_l(y))}{\partial x}B^y_l(y))+(B^x_i(x)\frac{\partial (B^x_i(x)B^y_j(y))}{\partial y},B^x_k(x)\frac{\partial (B^x_k(x)B^y_l(y))}{\partial y})\right) \) co dzięki strukturze produktu tensorowego funkcji B-spline, oraz dzięki temu że pochodna w kierunku y z funkcji B-spline zmiennej x daje 0, oraz pochodna w kierunku x z funkcji B-spline zmiennej y daje, dostajemy więc \( \sum_{i=1,...,N_x;j=1,...,N_y} dt*u^{t}_{i,j } \left( (\frac{\partial B^x_i(x)}{\partial x}B^y_j(y),\frac{\partial B^x_k(x)}{\partial x}B^y_l(y))+(B^x_i(x)\frac{\partial B^y_j(y)}{\partial y},B^x_k(x)\frac{\partial B^y_l(y)}{\partial y})\right) \).
3. Zmiany wywołane przez siłe działającą na system w czasie kroku czasowego \( (f_t,w)= (f_t(x,y),B^x_k(x)B^y_l(y)) \)

W metodzie explicite macierz układu to macierz masowa
\( \{ M_x \otimes M_y \}_{i,j,k,l}= \int B^x_i(x)B^x_k(x)B^y_j(y)B^y_l(y) = \int B^x_i(x)B^y_j(y)B^x_k(x)B^y_l(y) = {\bf M}_{i,j,k,l } \)
która posiada strukturę produktu Kroneckera i może ona zostać sfaktoryzowana w czasie liniowym przez solwer zmiennokierunkowy.
Wadą metody explicite jest fakt, iż pojedynczy krok czasowy nie może być za długi, inaczej metoda będzie niestabilna. Objawiać się to będzie eksplozją symulacji numerycznej. Wynika to z tak zwanego warunku Courant'a–Friedrichs'a–Lewy'ego (CFL condition). W ogólnym przypadku mówi on że symulacja numeryczna w której modelujemy pole skalarne \( u(x,t) \) jest stabilna jeśli \( C= \frac{u * dt}{ h} < C_{max } \) gdzie \( u = max u(x,t) \) oznacza maksymalną wartość symulowanego pola, \( dt \) oznacza długość kroku czasowego \( h \) oznacza rozmiar elementu siatki (co w przypadku analizy izogeometrycznej wynika z wartości punkt w wektorze węzłów), oraz \( C_{max} \) jest pewną stałą wynikającą z rozwiązywanego problemu. Praktyczne znaczenie tego warunku jest następujące. Jeśli będziemy zwiększać liczbę elementów siatki (zwiększać liczbę punktów w wektorze węzłów) wówczas \( h \) będzie się zmniejszać i żeby nasza symulacja nie "eksplodowała" musimy zmniejszać krok czasowy. Jest to duże ograniczenie, ponieważ oznacza to że zwiększanie dokładności w przestrzeni wymusza zmniejszanie rozmiarów kroków czasowych, co oznacza większy koszt symulacji. Pewnym rozwiązaniem tego problemu jest stosowanie symulacji implicite ze schematami umożliwiającymi stosowanie solwerów zmiennokierunkowych.